library(ggplot2)
library(dplyr)															######chromHMM染色质状态的p1
setwd("E:/5hmc_file/脑组织和GM12878的chromHMM数据")
anno=read.table("group意义.txt",head=T,sep="\t")
anno=data.frame(HMM.group=anno$STATE.NO.,description=anno$DESCRIPTION,num=anno$NUM,abbr=anno$MNEMONIC)

file=read.table("result_HMM_enrichment.csv",head=T,sep=",")
groupf=unique(file$con.file.name)

file1=file[file$con.file.name==groupf[2],]
file1$HMM.group=paste0("NO.",file1$HMM.group)
filee=merge(file1,anno,by="HMM.group")

col_names=c("HMM.group","OR","lower","upper","p.value")
result=data.frame(matrix(NA,15,ncol=5))
names(result)=col_names
group.HMM=unique(filee$HMM.group)
library(meta)
for(i in 1:15){
tdata=data.frame(HMM.group=filee$HMM.group,filee[,4:7])
tdata1=tdata[tdata$HMM.group==group.HMM[i],]
 tdata1$case_not=tdata1$case_in_region+tdata1$case_not
 tdata1$con_not=tdata1$con_in_region+tdata1$con_not
metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=tdata1,sm="OR",studlab = HMM.group) 
OR=exp(metaor3$TE.fixed)
upper=exp(metaor3$upper.fixed)
lower=exp(metaor3$lower.fixed)
result[i,]$HMM.group =group.HMM[i]
result[i,]$OR=exp(metaor3$TE.fixed)
result[i,]$upper=exp(metaor3$upper.fixed)
result[i,]$lower=exp(metaor3$lower.fixed)
result[i,]$p.value=metaor3$pval.fixed
}
result=merge(result,anno,by="HMM.group")
 result$FDR=p.adjust(result$p.value,method = "bonferroni")
 result$group="not significant enrichment"
 result[result$FDR<0.05&result$OR>1.5,]$group="significant enrichment"
 result$OR=log2(result$OR)
 result$lower=log2(result$lower)
 result$upper=log2(result$upper)
 result=result[order(result$num),]
 result1=result[1:9,]
 result2=result[10:15,]
 result=rbind(arrange(result1,group,OR),arrange(result2,group,OR))
 result$num=1:15
 result1=result

p1=ggplot(data=result,aes(y=OR,x=num,group=group,color=group))+geom_errorbar(aes(ymin=lower,ymax=upper),width=0.1,color="black")+
    geom_point(aes(shape=group),size=3)+coord_flip()+geom_hline(yintercept = log2(1.5),color="red",linetype="dashed")+scale_x_continuous(breaks = result$num,labels = result$description)+
    theme_light()+ylab("log2(OR)")+geom_hline(yintercept = 0,color="black")+
    scale_color_manual(values = alpha(c("#3C5488","#E64B35"),0.9))+theme(panel.grid.minor = element_blank(),legend.position = "none")	   
	   
	   
	 
													################组织特异性的图p2
													#其余画图前处理过程参见figure2文章用.R
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201204")
ASD=read.csv("meta.p.TWAS.ASD.statis.csv",header=T)	#这是阿尔兹海默症的数据
BD=read.csv("meta.p.TWAS.BD.statis.csv",header=T)
Depr=read.csv("meta.p.TWAS.Depr.statis.csv",header=T)
SCZ.BD=read.csv("meta.p.TWAS.SCZ.BD.statis.csv",header=T)
Depr=Depr[,-6]
ASD$group.source="Alzheimer"
BD$group.source="BD"
Depr$group.source="Depr"
SCZ.BD$group.source="scz.BD"
rt=rbind(ASD[1,],BD[1,],Depr[1,],SCZ.BD[1,])
rt$logpvalue=-log10(rt$P.value)
rt$counts=c(296,645,454,2417)
rt$group=ifelse(rt$P.value<0.05,"sig","nosig")
p2=ggplot(rt,aes(y=OR,x=logpvalue,color=group))+geom_point(size=3)+scale_color_manual(values = alpha(c('#3C5488','#E64B35'),0.8))+
    labs(y="OR",x="-log10 (P-value)",title="tis sepc enrichment")+
    theme_classic(base_size = 12)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+geom_hline(yintercept = 1.5,color="red",linetype="dashed")+theme(legend.position = "none")+
    geom_text_repel(data=subset(rt,rt$p.value<0.05),aes(label=tissname),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)

													###################TF富集图p4
setwd("E:/5hmc_file/motifbreakR_predict_enrich")
file=read.csv("enrichment_by_TF_predit_by_MBR.csv",head=T)
file=file[file$con.file=="nobias_AShM",]
file[file$OR=="Inf",]$OR=27
library(ggplot2)
library(ggrepel)
file$group="nosig"
file[file$FDR<0.1,]$group="sig"
file$logFDR=-log10(file$FDR)
p4=ggplot(file,aes(y=OR,x=logFDR,color=group))+geom_point(size=3)+scale_color_manual(values = alpha(c('#3C5488','#E64B35'),0.7))+
    labs(y="OR",x="-log10 (FDR)",title="TF enrichment")+
    theme_classic(base_size = 12)+geom_vline(xintercept = -log10(0.1),color="red",linetype="dashed")+geom_hline(yintercept = 1.5,color="red",linetype="dashed")+theme(legend.position = "none")+
    geom_text_repel(data=subset(file,file$FDR<0.1),aes(label=TFname),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)


  
#############################################################调控特征出现的可能性和羟甲基化升降一致性的统计
setwd("E:/5hmc_file/H3k的分析")
sel=list.files(pattern = "Brain")
library(ggplot2)
library(ggrepel)
library(RColorBrewer)

alias=c("AG","AC","GM","HM","ITL","MFL","SN","FBF","FBM")
for(i in 2:length(sel)){
ftp=read.csv(sel[i],head=T)
names(ftp)=c("term",names(ftp)[-c(1)])
ftp$group=alias[i]
file=rbind(file,ftp)
}
file$logpvalue=-log10(file$pvalue)

getPalette = colorRampPalette(brewer.pal(9, "Set1"))

p5=ggplot(file,aes(y=same_ratio,x=logpvalue,color=term))+geom_point(size=3)+scale_color_manual(values = alpha(c('#C71585','#D2691E','#E64B35','#3CB371','#4169E1','#708090','#000000'),0.9))+
  labs(y="same_ratio",x="-log10 (Pvalue)",title="H3k")+
  theme_classic(base_size = 12)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
  geom_text_repel(data=subset(file,file$pvalue<0.05),aes(label=group),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)

########################
########################
layout <- matrix(c(1, 2, 3, 4), 2,2,byrow = T)
multiplot(plotlist=list(p1,p2,p5,p4),layout=layout)